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We study a discrete nonlinear Schrodinger lattice with a parabolic trapping potential. The model, describing, 
e.g., an array of repulsive Bose-Einstein condensate droplets confined in the wells of an optical lattice, is ana- 
lytically and numerically investigated. Starting from the linear limit of the problem, we use global bifurcation 
theory to rigorously prove that - in the discrete regime - all linear states lead to nonlinear generalizations 
thereof, which assume the form of a chain of discrete dark solitons (as the density increases). The stability of 
the ensuing nonlinear states is studied and it is found that the ground state is stable, while the excited states 
feature a chain of stability/instability bands. We illustrate the mechanisms under which discreteness desta- 
bilizes the dark-soliton configurations, which become stable only inside the continuum regime. Continuation 
from the anti-continuum limit is also considered, and a rich bifurcation structure is revealed. 



I. INTRODUCTION 

The experimental realization of atomic Bose-Einstein 
condensates (BECs)^ has triggered an intense activity in 
the study of purely quantum systems at almost macro- 
scopic scales. From a theoretical standpoint, many effects 
related to BEC physics can be described by lowest-order 
mean-field theory, namely the Gross-Pitaevskii equation 
(GPE)^!!. The latter is a nonlinear Schrodinger (NLS) 
equation, incorporating an external trapping potential, 
with the nonlinearity effectively accounting for inter- 
atomic interactions. In the absence of the nonlinear 
term, the GPE becomes a linear Schrodinger equation 
for a confined single-particle state; in this limit, and in 
the case of, e.g., a harmonic external potential, the lin- 
ear problem becomes the equation for the quantum har- 
monic oscillator characterized by discrete energies and 
corresponding eigenstates^. Then, one may generalize 
this simple physical picture, and regard the GPE as a 
model for a self-interacting macroscopic quantum oscilla- 
tor, in such a case, the use of analytical and/or numerical 
techniques for the continuation of these linear eigenstates 
(supported by the particular type of the external trap- 
ping potential) leads to purely nonlinear states of the self- 
interacting quantum oscillator. Such nonlinear states can 
be found both in one-dimensional (ID)^"^ and higher- 
dimensional settingsSrJ^. Notice that in the ID setting, 
and for BECs with repulsive interatomic interactions, the 
nonlinear states assume the form of dark solitons, which 
have been studied extensively both in nonlinear optic&i^ 
and the physics of atomic BECs^. 

In this work, we consider and analyze the discrete ver- 
sion of the GPE model, namely a discrete NLS (DNLS) 
equation}^ , which incorporates a (discrete) harmonic po- 
tential. This model is motivated by the physical setting 
of a BEC confined in a combined non-negligible harmonic 
trap and strong periodic potential, the so-called optical 



lattice, where rich physical properties and nonlinear dy- 
namics have been revealed^^"— . Optical lattices are 
generated by a pair of laser beams forming a standing 
wave which induces a periodic potential; thus, for a BEC 
loaded in an optical lattice, the trapping potential in the 
GPE can be regarded as a superposition of a harmonic 
trap and a periodic potential. If the harmonic potential 
is very weak as compared to the optical lattice, it can ap- 
proximately be ignored; then, the stationary states of the 
GPE (which includes solely the periodic potential) can be 
found in the form of nonlinear Bloch waves, which have 
the periodicity of the optical lattice (see, e.g., Ch. 6 in 
Ref.-^ and references therein). In the same case (i.e., in 
the absence of the harmonic potential) , if the optical lat- 
tice is suHiciently deep (compared to the chemical poten- 
tial), the strongly spatially localized wavefunctions at the 
lattice sites can be approximated by Wannier functions 
(see, e.g., Ref.^") and the tight-binding approximation 
can be applied; then, the continuous GPE is reduced to 
the DNLS equation^ '^^i^^'^" , a model which has already a 
long history in the physics and mathematics of nonlinear 
lattices^^. Notice that the validity of this model assumes 
intra-well phase coherence and, thus, it cannot be used in 
situations such as the superfluid-to-Mott insulator phase 
transitional or, generally, when strong correlation effects 
come into play (see, e.g., the review^). Nevertheless, the 
model under consideration, apart from being motivated 
by the physics of BECs loaded in optical lattices - where 
it can be regarded as a macroscopic quantum harmonic 
oscillator on a lattice - it may also apply in other physical 
settings, including discrete nonlinear optics^ and nonlin- 
ear lattice theoriesi^. 

Here, our scope is to study the existence, bifurcations 
and stability of nonlinear states emerging in this setting 
for values of the lattice spacing a ranging from the dis- 
crete regime (a — 0(1)) to the so-called anti-continuum 
(AC) limit (a — oo). First, we revisit the linear limit of 
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the problem (studied some time ago in Ref.^*) and gen- 
eralize the corresponding linear considerations towards 
finding analytically and numerically the nonlinear states 
supported by the system. Then, we use global bifur- 
cation theory Sl ii H to rigorously prove that - in the dis- 
crete regime - each linear eigenstate of the system can 
be continued to a nonlinear counterpart. This way, we 
find all such nonlinear states, namely the ground state, as 
well as excited states which, within the strongly nonlinear 
regime, acquire the form of a chain of dark solitons. We 
also study the continuation from the AC limit, through 
a detailed numerical bifurcation analysis, and find that 
there exist states without a linear counterpart. 

Furthermore, we study the stability of the nonlinear 
states in the framework of the linear stability analysis 
[so-called, Bogoliubov - de Gennes (BdG) analysis, in 
the realm of BECs] focusing on the effect of discreteness. 
We reveal a fundamental difference of the discrete sys- 
tem and its continuum counterpart: we find that the dis- 
creteness renders the excited states more unstable (only 
the ground state is found to be always stable) through 
the emergence of instability bands; the pertinent band 
structure depends on the lattice spacing and the non- 
linearity strength (as measured by the chemical poten- 
tial fi). Contrary, in the continuum case, not only the 
ground state but even some of the excited states (such 
as the first and second ones) are stable deeply inside the 
nonlinear regime^^. We also perform numerical simula- 
tions to follow the evolution of the first two (unstable) 
excited states, namely of the single discrete dark soliton 
and of the dark soliton pair. We show that the (oscil- 
latory) instability thereof manifests itself by setting a 
quiescent dark soliton configuration into an oscillatory 
motion, which can be explained by resonance effects be- 
tween the eigenfrequencies of the soliton modes and the 
intrinsic excitation frequencies of the underlying system. 

The paper is organized as follows. In Section II, we 
present the model as motivated by the physics of BECs 
loaded in optical lattices, although, as indicated above, 
our considerations can be relevant to other fields of ap- 
plications such as nonlinear optics. In Section III, we 
also study analytically and numerically the linear limit 
of the model. In Section IV, we consider the fully non- 
linear problem and show, in particular, how nonlinear 
eigenstates emerge from linear ones; the anti-continuum 
limit of the system is also studied. In Section V, we ana- 
lyze the stability of the excited nonlinear states (i.e., the 
single dark soliton and the two dark soliton states) and, 
finally, in Section VI, we summarize our conclusions. 



II. MODEL AND METHODS 

A. Physical motivation and the model 

We consider an atomic BEC confined in a highly 
anisotropic harmonic trap, Vht-, with frequencies ujx and 
u;± = ujy = Wz, such that lOx ^ In the mean- 



field approximation, and for sufficiently low temperatures 
(so that thermal and quantum fiuctuations can be ne- 
glected), the BEC dynamics can be described by the fol- 
lowing effectively one-dimensional (ID) GPE^, 

where ^'(a;,t) is the macroscopic BEC wavefunction nor- 
malized to the number of atoms, namely J j^'pda; = N, 
while (71D — 2hujj^a is the effectively ID coupling con- 
stant, with m being the atomic mass and a the s-wave 
scattering length, assumed to be positive (i.e., the inter- 
atomic interactions are repulsive). Finally, the external 
potential, Vext (x) , in Eq. ([T]) takes the form 

Vext{x) = VHT{x) ^^mujlx^. (2) 

Equation ([1]) can be expressed in the following dimen- 
sionless form: 

where — 2a|\E'p is the normalized density, while 
length, time and energy are respectively measured in 
units of 2a, aj_, wj^ and hijjj_; the potential V{x) in 
Eq. ([3]) is given by: 

Vix) ^ ^n'x', (4) 

where Q = u;x/u!± is the normalized harmonic trap 
strength. Notice that apart from the BEC context, 
Eq. appears also in studies in the nonlinear optics con- 
text (see, e.g., Ref;^): there, ^' is the normalized electric 
field envelope, t denotes the propagation direction, while 
the parameter fl accounts for the change in the refractive 
index of the medium in the x-direction (transverse to the 
propagation). 

In our analysis below, we consider the discretized ver- 
sion of Eq. ([3]), namely the following DNLS model: 

= -^^2^, + if^'(«jf + IV'jf V-,, (5) 

where the overdot denotes time derivative, a is the lattice 
spacing, and A2^pj = ipj+i — 2'0j + tpj-i is the discrete 
Laplacian. In the continuum limit of a — > 0, Eq. ([5]) 
is reduced to the continuum CP model of Eq. ([3]). No- 
tice that in the absence of the nonlinear term, Eq. ([5]) 
is the time-dependent problem for a quantum harmonic 
oscillator (QHO) on a lattice. On the other hand, in the 
presence of the nonlinear term, the model can be consid- 
ered as a self-interacting macroscopic QHO on a lattice. 

The DNLS of Eq. ([5]), apart from being interesting in 
its own right, is also relevant to the physics of atomic 
BECs confined in strong optical lattices. To further 
elaborate on the above, let us assume that the exter- 
nal potential V^xt in Eq- © incorporates a periodic 
optical lattice potential, Vqli created by two counter- 
propagating laser beams of wavelength A^^'^^ ; in such a 
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case, Vext{x) = Vht{x) + Vol{x), with the optical lattice 
potential being given by: 



Vol{x) = Vosin^{kx), 



(6) 



where Vq and k = 2n/X denote the strength and 
wavenumber of the optical lattice, respectively. Then, 
considering the special case of a sufficiently strong opti- 
cal lattice, such that Vb ^ /x (where fi is the chemical 
potential), we may follow the analysis of Ref.^ and show 
that Eq. (H)) can be approximated by a DNLS model for 
the wavefunctions V-'j (t) in the different wells (denoted by 
the index j). Particularly, we assume that the effective 
harmonic trap frequency at each well, uJx = \j2V()^^^Tn, 
is such that Cox ^ oJx, we may employ the tight-binding 
approximation and decompose the BEC wavefunction 
^{x, t) as a sum of the wavefunctions <^j{x — Xj) localized 
around the center of each well, namely. 



m{x,t) = Y,^,{t)^,{, 



(7) 



where the wavefunctions are normalized to unity, and 
the total number of atoms in the condensate now reads 
TV = "Y^- Nj — J2 IV-'j P (where Nj is the number of atoms 
at the well j). Substituting Eq. © into Eq. (P), muhi- 
plying by and integrating over x, yields the following 
equation for the wavefunctions ipj (see Refsi^i^), 



(8) 



To derive the above equation, we have used the (quasi) 



orthogonality relation j dx^i^j 



we have kept 



only terms including spatial integrals of first-neighbor 
wavefunctions, and we neglected terms proportional to 
/ dx^'j^^j^ and / dx^^^j±i. The constants K and Ej 
in Eq. Q are given by: 



K ^ 



dx 



dx 



2m dx dx 



'^jVext{x)<^j±l 



2m 



dx 



(9) 
(10) 



while g — ^id / Equation ([5]) can readily be 

made dimensionless measuring length, time and energy 
in units of the lattice spacing a = A/2, = H/El, and 
El = 2Efi = fi? /ma^ (where Eji is the recoil energy), 
respectively. In these units, and employing the transfor- 
mation, 



ibj — -1 / — exp 
V 5 



nujT, 



(11) 



where e = K/ El and VL = lJx/lol^ respectively. 

Formally speaking, the DNLS Eq. (fT2|) is a variant 
of Eq. ([SJ, but there are also some differences arising 
from the dependence of the trap strengths and coeffi- 
cients of the kinetic terms on the lattice spacing a. From 
a physical viewpoint, Eq. (|12|) applies for the regime 
corresponding to moderate values of a: this is due the 
fact that the assumptions for the derivation of Eq. (fT2|) 
become invalid for small or large values of the lattice 
spacing. Nevertheless, it can be found that there exists 
a certain range of a-values (for a given harmonic trap 
strength fi), where Eq. is equivalent to Eq. ([SJ - 
the formal discretization of Eq. using experimen- 
tally relevant parameters^S, for a rubidium condensate 
confined in a trap with frequencies = 27r x 10 Hz and 
= 27r X 100 Hz, and total number of atoms N k, 2000, 
one may find that the ratio Cl'^/K takes values in the in- 
terval 0.001 < n'^/K < 0.01. Accordingly, for the fixed 
value of the normalized trap strength = 0.1 (which will 
be used below), if the lattice spacing takes values in the 
interval 0.25 < a < 0.6, then Eqs. (O and ([H]) become 
equivalent. 

Thus, the model Eq. ([S]) is related to the continuum 
model Eq. ([3]) (for a — > 0), describing dynamics of 
harmonically confined BECs or dynamics of beams in 
graded-index waveguides, while it can also be used - in 
the strongly discrete regime (a < 1) - to describe the 
dynamics of arrays of BECs in optical lattices. 

It should be noted in passing that, in what follows in 
our analysis, as the number of atoms tends to zero, quan- 
tum effects considered in a number of recent works^i"— 
become important; in such a case, applicability of the 
mean-field approximation becomes questionable. Never- 
theless, our aim here is to utilize the model at hand - as 
a relevant mathematical limit ~ which can be explored 
to identify the nonlinear states emerging from the lin- 
ear ones in the regime where the mean-field description 
is the appropriate one (i.e., for sufficiently large atom 
numbers). 



B. Stability analysis approach 

Below, we will present results concerning the stability 
of nonlinear states of Eq. ([5|). In fact, we will perform 
a linear stability analysis based on the so-called (in the 
context of BECs) BdG equations^'-. In particular, once 
a real, stationary state, is identified by means of a 
fixed point algorithm (e.g., a Newton- Raphson method), 
we consider small perturbations of this state of the form, 

Mt) = l^f + (".e-^"* + ^^*e'"*)]e-'^^ (13) 



we can express Eq. ([S]) as follows: 



(12) 



where the asterisk denotes complex conjugation. Substi- 
tuting this ansatz into Eq. ([S]), and linearizing with re- 
spect to Uj and Vj, we obtain the linear stability (BdG) 
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equations 



i?-/i + 2|V' 



(0)|2 



(0)|2 



-UJVj 



(14) 
(15) 



where H — — (l/2a^)A2 + iri^(aj)^ is the single-particle 
operator. Solving these equations one can find the eigen- 
frequencies uj = ujr + iuji and the amplitudes Uj and Vj 
of the normal modes of the system. Note that due to 
the Hamiltonian nature of the system, if uj is an eigenfre- 
quency of the Bogoliubov spectrum, so are — w, uj* and 
—UJ*. A stable (unstable) configuration corresponds to 
io, = {uj, ^ 0). 

An important quantity resulting from the BdG analy- 
sis is the energy carried by the normal mode with eigen- 
frequency uj. This is given by the following expression: 



E 



dxi\u\^ - \v\'^)uj. 



(16) 



In brief, the energy measure of (1161) yields the energy 
difference between a perturbed state and an equilibrium 
(fixed point state), as an explicit calculation (see equa- 
tions (5.73)-(5.77) of the Ref.^) clearly illustrates. The 
sign of this quantity, known as Krein sigv?^-, is a topo- 
logical property. Importantly, if the normal mode eigen- 
frequencies with opposite energy (Krein) signs are in res- 
onance then, typically, there appear complex frequencies 
in the excitation spectrum, i.e., a dynamical instability 
occurs^i. In order to further elaborate on such a pos- 
sibility, we note that modes with complex or imaginary 
frequencies carry zero energy, while anomalous modes - 
associated with the presence of dark solitons in the con- 
figuration - have negative energy (see, e.g., Sec. 5.6 of 
Ref.^). The presence of anomalous modes in the excita- 
tion spectrum is a direct signature of an energetic insta- 
bility or, in other words, is an evidence that the station- 
ary state over which the BdG analysis is applied is not 
the ground state of the system. 

The above analysis scheme will be used in Sees. IV. B 
and V below. 



III. THE LINEAR PROBLEM 

Let us start our analysis by considering at first the 
linear counterpart of Eq. ([5]) resulting from the substitu- 
tion t/jj — )• ipj exp{—iEt) (where E denotes the energy), 
namely. 



1 



1 



—A2ijj + -n'iajy^,=E^,. 



(17) 



The above equation is the discrete version of the eigen- 
value problem describing a QHO on a lattice. The en- 
ergy spectrum, as well as the profiles of the pertinent 
eigenstates, will be used below to construct (numerically) 
solutions of the full nonlinear problem. Following the 
methodology devised in Ref.^^, it is possible to solve the 



discrete QHO eigenvalue problem by considering the fol- 
lowing (continuous) Hamiltonian operator acting on the 
wavefunction ^'(a;,t): 



^n^x^^ = E'^, (18) 



where p = —id/dx is the momentum operator, A is the 
tight-binding constant, a is a constant, and E' is the 
energy. Equation (fT^ is identical to Eq. (IT71) in a dis- 
crete coordinate space: indeed, letting \I'(a;,i) ^j{t)i 
X — >■ aj and V[x) — >■ ^f2^(aj)^, the translation operators 
exp(±jQ;p) act on the wavefunctions as ex.'p{±iap)^ j = 
il>j±i and a corresponds to the lattice spacing. Then, 
adding on both sides of Eq. (IT8|) the term 2Aij)j , we find 



Mi'j+i + V'j-i 



2^,j) + Ul\ajf^,^Eij^, (19) 



where we need to identify A = l/2a^ and E ^ E' + 2A. 
We thus need to solve the continuous eigenvalue problem 
of Eq. p^ . This can be done by expressing it in the 
momentum representation (i.e., p = p and x = id/dp), 
where it can be written as the following Mathieu-type 
equatio: 



,24 



+ [b-2qcos{2v)]<j>{v) ^0, (20) 



ap 8A , 8E' 



where 



and (f){v) is the Fourier transform of Equation ((20)) 
possesses a well-known energy spectrum and solutions 
(see, e.g., Ref.— ). Projecting these solutions on the 
Hilbert space where x/a = n {n = 1,2, ...) will provide us 
with the solutions of Eq. The solutions of Eq. ([20]), 

namely, 



^(f-»)(„) = iV,ce2„(t;;g), 



n = 0,l,2,--. (22) 
n = 1,2,3,- •• (23) 



are called Mathieu functions and arc periodic, of period 
TT. For different values of q, these solutions correspond 
to characteristic values of b, namely A2n and B2n for the 
even and odd eigenfunctions, respectively, from which we 
deduce the following energy spectrum: 



(24) 
(25) 



for the even and odd eigenfunctions, respectively. The 
energy spectrum has a simple form, both in the con- 
tinuum limit, corresponding to a — >■ 0, and the anti- 
continuum limit, corresponding to a — ?■ oo; the respec- 
tive analytical expressions for A2n{Q) and B2niQ) can be 
found in Ref."^^. As is expected, in the continuum limit. 
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FIG. 1. (Color online) The ground-state energy Eo for the 
linear quantum harmonic oscillator as a function of the lat- 
tice spacing a, for two values of the trap frequency, f2 = 0.1 
and Q = 0.05. Solid and dotted lines show the energy spec- 
trum as found by solving the QHO eigenvalue problem, while 
(red) circles and (blue) crosses show the respective solutions 
obtained from the Mathieu equation (|20|l . The dashed and 
dashed-dotted lines show the analytical result of Eq. (|28|l . 
The inset shows the wavefunction profile for a — 10. 

we recover the familiar equidistant QHO energy spectrum 
virith energies 

E:^=(n+^^n, n = 0,l,2,..-, (26) 

while in the anti-continuum limit, the energy spectrum 
becomes parabolic and has the form 

E'^^^^n^a^n^, n = 0,l,2,.... (27) 

The eigenfunctions in the coordinate space can be ob- 
tained upon Fourier transforming the Mathieu functions 
of Eqs. ([22])- ([231). As shown in Ref.— , these eigenfunc- 
tions are very similar to the Hermite polynomials and 
coincide with the latter in the continimm limit. 

The analytical results presented above can directly be 
compared with numerics. We first study the ground state 
energy, for two different oscillator frequencies, spanning 
all the allowable range of values of a, from the continuum 
limit (a — ^ 0) to the anti-continuum one (a — > oo). The 
energy given by Eq. (I24p has an approximate analytical 
form, namely, 

for sufficiently small values of q. 

In Fig. [1] we compare the dependence of the ground- 
state energy on the lattice spacing a found by numeri- 
cally solving the QHO eigenvalue problem, with the one 
obtained by solving the Mathieu equation ([20]) : we also 
show the approximate analytical result of Eq. (IMl) . The 
analytical solution is only a good approximation for suf- 
ficiently large a - or for small values of the parameter 
q. As expected, the ground state energy in the contin- 
uum limit, a — > 0, is equal to 51/2, while in the anti- 
continuum limit, a oo, the energy is independent of 
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FIG. 2. (Color online) Top right panel: Energy of the lowest 
four excited states as functions of the lattice spacing a, for a 
trap strength Q — 0.1. Solid lines (blue, red, green, magenta, 
correspond to the 1st-, 2nd-, 3d-, and 4th-excited states, re- 
spectively) indicate the energy obtained by solving the QHO 
eigenvalue problem, dashed lines show the respective solutions 
obtained from the Mathieu equation (|20p . and circles show 
the respective analytical results of Eq. (|25p . Top left panels: 
spatial profiles of the 1st-, 2nd-, 3d-, and 4th-excited states 
for Q = 10 (corresponding to the anti-continuum limit); solid 
lines depict the 1st- (top) and 3d- (bottom) excited states, 
while dotted lines depict the 2nd- (top) and 4th- (bottom) ex- 
cited states, respectively. Bottom panels (from left to right): 
spatial profiles of the 1st-, 2nd-, 3d-, and 4th-excited states 
for a — 1 (corresponding to the discrete regime). 



the trap strength [cf. Eq. (P7)) ]. The latter result can 
be understood from the profile of the wave function in the 
anti-continuum limit, as seen in Fig.[T] In this limit, the 
only excited site is j = which, according to Eq. (|17l) . 
yields E = and asymptotically goes to zero. 

Next, we study the spectrum of the first four excited 
states. The equidistant spectrum in the continuum limit 
- see Fig. [5]- becomes parabolic in the fully discrete case. 
In the anti-continuum limit it is observed that the ex- 
cited states become degenerate in pairs. Again the profile 
of the wave functions explains this result: in this limit, 
e — 0, as seen from Eq. (|17l) the energy depends solely on 
the potential which is quadratic. As shown in the inset in 
the top panel of Fig. [5] (where the profiles of the first four 
excited states are shown for a = 10), the energy needed 
to excite symmetrically or anti-symmetrically (with re- 
spect to the center) the first neighboring sites is exactly 
the same due to the quadratic nature of the potential. 
On the other hand, in the discrete regime, the wavefunc- 
tion profiles are characterized by the number of nodes, 
i.e., n- nodes for the n-th excited state; pertinent profiles, 
for the first four excited states, are shown in the bottom 
panels of Fig. [5] (for a = 1, and the same trap strength, 
fl = 0.1). 
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IV. EXISTENCE AND BIFURCATIONS OF SOLUTIONS 
IN THE FULLY NONLINEAR PROBLEM 

A. Continuation from the linear to the nonlinear regime 

In this section we will study the fully nonlinear case. 
Our analysis considers an arbitrary number of K + 2 os- 
cillators equidistantly occupying an interval [—L, L] of 
length 2L, with spacing a — -jf^- Thus, the oscil- 
lators are occupying the points Xj — —L + ja, j = 
0,1,2, K + 1 of the interval [— L, L], discretized as 



L = xo < a;i < . . . < xk+i = L. 



(29) 



We consider the case of real discrete wavefunctions. 
For the discrete wavefunctions at each point Xj, j = 
0, . . . , K + 1, of (P^ . we use the standard shorthand 
notation ip{xj) = In some cases we shall also use 
the shorthand notation ip for the vectors of M^+^, i.e., 

^ ■■= mf=V- 

First, we use the transformation -tpj — >■ ipj exp(— i/it) 
(where fi is the chemical potential) to reduce Eq. ([5]) to 
its time-independent counterpart, 

- ^A2Vj + lnHaj)'^P, + I^.IV, = mV'j, (30) 

for j = 1, . . . , K, and assume that the wavefunctions are 
satisfying Dirichlet boundary conditions at the endpoints 
xq = —L and xk+i — L, namely: 



K+l 



= 0. 



(31) 



Our aim is to use the solutions of the linear problem 
obtained in the previous section in order to find pertinent 
solutions in the nonlinear regime. Our analysis starts 
by proving that the energy spectrum of the nonlinear 
equation arises from the relevant spectrum found in the 
linear case. Before proceeding further, it is relevant to 
note that Eq. pOI) . with the boundary conditions (I5T]) . 
possesses two conserved quantities: the Hamiltonian H 
(the energy of the system) and the number of atoms N, 
respectively given by: 



1 ^ 



K 



(32) 
(33) 



Notice that the presence of the prefactor a in the defi- 
nitions of H and N suggests that in the limit of a — > 
Eqs. (|32|) and (pS)) provide the respective Hamiltonian 
and number of atoms of the continuum GPE, Eq. ([3]). 

The continuation to the nonlinear regime from the 
linear states (IT7| . i.e., the bifurcations of solutions of 
Eq. (1301) from solutions of the corresponding linear prob- 
lem, cf. Eq. P7)) . can be justified analytically by using 
global bifurcation theory - see Refs.-^ and Section 15.7 



of Ref.^^. In this setting, we will apply the global bifur- 
cation theorem of Rabinowitz (see Theorem 1.3, p. 490, 
of Ref.^ and Theorem 15. C, p. 668, of Ref.-^^), which we 
now recall for reasons of completeness 

Theorem 1 Assume that X is a Banach space with 
norm \ \ ■ \ \x- Consider the map J-{p., ■) : X ^ X, S M, 
J^(/i, •) = /i£ • +'H(/i, •), where C : X ~> X is a compact 
linear map and 'H{^, ■) : X ^ X is compact and satisfies 



lim 

||ti||x^O 



||H(/i,w)|| 



X 



0. 



\u\\x 



(34) 



// is a simple eigenvalue of C, then the closure of the 
set 

C = {{n,u) e M X X : (/i,u) 
solves u — J^(/i, u) = Q, 0}, 

possesses a maximal continuum (i.e. connected branch) 
of solutions C which branches out of {X* , 0) and C either: 

(i) meets infinity m M x X or, 

(ii) meets u — in a point (/i,0) where p, ^ X* and 4 
is an eigenvalue of C. 

To apply Theorem [U we need some preparations, in or- 
der to rewrite ([30)) in the form tp — ^C{tp) + 'H{n,^) =0 
requested by the theorem. As a first step, we will de- 
fine and discuss the properties of the linear operator C 
through the linear eigenvalue problem (IT7l) - ([3T|) . which is 
the eigenvalue problem for the linear operator 



1 

2a 



(35) 



J = 1, . . . , K , supplemented with the Dirichlet boundary 
conditions (I5T]) . We shall also discuss some properties of 
the eigensolutions of (135]) related to the number of sign- 
changes of the discrete eigenfunctions. These properties 
will be useful for distinguishing between the possibilities 
(i) and (ii) described by Theorem [T] As a second step, 
we will define the nonlinear operator T-L through the non- 
linearity of (PH)) . 

The operator psp is strongly positive and selfadjoint 
on the Hilbert space 



j=K+l TaK+2 



V'O = i'K+l = 0}, 



having the role of the Banach space X which is referred 
in Theorem [TJ The Hilbert space X is endowed with the 
norm ([33]) . i.e.. 



The operator K possesses K simple eigenvalues, 

Q<Eo<E2< ...<Ek-i. (36) 

Furthermore, the Krein-Rutman theorem (see p. 122 of 
Ref,^, and Section 7.8, p. 289 of Ref,^), implies that 
the principal eigenstate associated to the principal eigen- 
value Eq is positive in the sense that > Q for all 
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j = 0, . . . , iiT + 1, '0 has at least one positive coordinate 
and satisfies the boundary conditions ([5Tt . On the other 
hand, the eigenvalue problem for the operator (|35|) with 
the boundary conditions ([5T|) . is the discrete analogue of 
the Sturm-Liouville problem for the QHO 

- + ^n^x^ilj{x) ^Xtp, -L<x<L (37) 

i^i-L) = iPiL) = 0, (38) 

for which the classical Sturm-Liouville theorem holds 
(see p. 454 of Ref.-^). For instance, (I37l)-(l38l) has a 
countable sequence of eigenvalues Ao < A2 < . . . , with 
corresponding eigenfunctions 'ipo{x),ipi{x), . . . , and the 
eigenfunction '0„(x), n = 0,1,..., has exactly n ze- 
ros (nodal points) on {—L,L). Continuing the discus- 
sion from the end of Section IIIIl the discrete eigenstates 
-0", n — 0,1, . . . , K — 1 corresponding to the eigenval- 
ues (j36p interpolate the continuous eigenfunctions uq(x), 
n — 0, 1, ... , K, and they have exactly n nodal points, 
n = 0,1, . . . , K — 1. We remark that for any e > the 
interpolation and the "nodal properties" of the discrete 
eigenfunctions of ((35)) have been described in RefSi ^^'^^ . 
Under this observation, for each n — 0, . . . , K — 1, we 
may define the following sets in X, 

Sn ■■= - {^j}jZo^' e K'^+' : 00 = V'K+i = 0, 
and has exactly n nodal points}. (39) 

The sets Sn are clearly open in X, since for any £ 
Sn we may construct an r-neighborhood B[i}),r) :— 
{4> £ X : 110 — 0||x < r}, lying in Sn, by considering r 
sufficiently small. For instance, for r-sufficiently small we 
get sufficiently small perturbations of the coordinates of 
in Sn and thus all the vectors of X being in B{ij), r), 
have the same number of nodal points (i.e., a small per- 
turbation oi £ Sn lies in Sn)- 

To conclude with our preparations, we write the non- 
linear eigenvalue problem (IBOp in the form of an operator 
equation in X, as follows; 

r(V') - mV' + •F(V') = 0, V e X, (40) 

where J- : X ^ X is the nonlinear operator 

J-(0;), = |V,|V,- 

The linear operator T : ^ ^ is invertible. Its inverse 
T^^ : X ^ X \s also symmetric and it readily follows 
that Vn '■= -g-, n — 0,1, . . . ,K —1, are also simple eigen- 
values of : X ~^ X. We may write Eq. (|40|) as 

^- ^JiT'^{1i^)+T-^F{^)^0 (41) 

Equation (I41[) is actually in the form requested by The- 
orem [T] 

0-M/:(V')+H(0) = o, (42) 

with the linear operator C := T^^ : X ^ X and the 
nonlinear operator H :— T~^T : X ^ X being compact 
since they are acting on the finite dimensional space X. 



The map F is defined by the cubic nonlinearity and, thus, 
it is not difficult to verify that % satisfies condition (|34|) 
of Theorem [TJ Hence, all the assumptions of Theorem 
[1] are satisfied, justifying that {En, 0), n = 0,1, . . . ,K — 
1 is a bifurcation point for the problem h3(f(l . We may 
summarize in the following: 

Proposition 1 For any a > 0, there exists a maximal 
continuum (i.e. connected branch) of solutions Ce„ of 
Eq. iSOp . n — 0,1, ... , K-^l, bifurcating from {En, 0) and 
Ce„ either (i) meets infinity inWxX , or (ii) meets = 
in a point {fi, 0) where fi 7^ En and i is an eigenvalue of 
£. 

We proceed by discussing some geometric properties 
of the branches Ce„- Considering the eigenstates -0" 
of Eq. p?)l corresponding to the eigenvalues £"„ , the lo- 
cal bifurcation theory and the implicit function theorem 
[see^e (Theorem 13.4 pg. 171 & Theorem 13.5, pg. 173)] 
guarantees that the local branch Ce„ can be locally rep- 
resented by the curve 

{^i,^p) : {-S,S) -^RxX, 

for some S sufficiently small, such that 

m=En, X(0) = 0, 

{fi{s),^{s))^{l,{s),s{r+x{s))), \s\<S, (43) 

where ||x(s)||x = 0(|s|), in the neighborhood of the bi- 
furcation point {En,0). Furthermore, there is a neigh- 
borhood of {En,0), such that any zero of the equation 
lies on this curve, or is of the form {En, 0)). 

Proposition 2 Consider the local representation ^S\ l of 
the branch Ce^- Then, /i'(0) = 0, /x"(0) > and the 
branch is locally concave up. 

Proof: We insert the expression ((^(s), ■0(s)) = 
{ijl{s) , s4'"' + sx{s)) in Eq. ([50]) and we divide by s. 
Then we obtain the equation (recalling that Xo(s) = 
XK+i{s) = 0), 

i^{^W; + xM) = -^^2(v; + xA^)) 

+ ]p^{ajf{^-+xAs)) 

+s>,"+X,(s)P(^7 + X,(s)). 

(44) 

We now differentiate Eq. pi)) with respect to s, namely, 
^i\s){^p-+xA^))+^,{s)^As) = -^A2xK^) 
+ ]p\ajfx'As) + 2.|^; + x,(^)l'(0; + XM) 

+3sV" + X,(«)Px;(^), (45) 

and by setting s = in Eq. (|^5]) and using Eq. P5)) . we 
have: 

-^A,x;(o) + ii^^(a,rx;(o) 

= /x'(0)V;+i?„x;(0). (46) 



8 



Multiplication of Eq. p6|) by -0" and summation by parts, 
yields 



K+l 



K+l 



j=0 J=0 
-, K+l K+l 

K+l K+l 

= ^ m'(0)|V;P + J2 x'{0)jE^^7 (47) 



J=0 i=0 

Since En and i/"? solve Eq. ([T7]), we have that 



K+l 



K+l 



2a 



j=0 j=0 

j=o 

Thus, from Eq. (gT]) we get that 

K+l 

implying that /i'(0) = 0. Next, by differentiating Eq. pS)) 
with respect to s, and setting s = 0, one obtains the 
equation 

-^Ax;'(0) + ^n^ajfx'^iO) + 2\r,\^r, 
Working as before, this time we derive that 

K+l K+l 

2Ei^"i' = ^"(o)Eiv^"i'- 

i=o i=o 

Therefore, ^"(0) > 0. o 

Proposition [21 actually states that at least locally, the 
graph of the function \i{N) iC^ curve) is concave up 
and monotone, thus locally invertible. By interchanging 
the axes and plotting the graph of as a function of 
/i, we recover that has locally the same concav- 

ity properties (see, e.g., Ch. 13 of Ref.^^ for bifurcation 
diagrams), as stated in Proposition [21 

It remains to show that the branches Ce^ are global, 
i.e., that the option (ii) of Theorem{l\ should be excluded. 

Theorem 2 For any a > 0, the maximal continuum 
(connected branch) of solutions Ce„ of Eq. I130\) bifurcat- 
ing from {En, 0) meets infinity in K. It is locally concave 
up and is not possessing a maximum (minimum) point. 

Proof: (a) Recall that any solution {fi, V') close to (£'„, 0) 
has the same number of nodal points as the eigenstate ■0" 
corresponding to the eigenvalue En. This is due to the 




FIG. 3. The number of atoms A'^ as a function of the chemical 
potential jj, (for a = 0.8 and = 0.1) for the three lowest 
states: the ground state (solid line), the first excited state 
(dashed line), and the second excited state (dotted line). Each 
branch begins from the linear limit (A'^ = 0), where [i equals 
the energy of the corresponding linear state. The insets show 
the profiles of these nonlinear states for [i, = 1.2. 



C^-representation of the solution -0 as 0j(s) — sip" + 
sxjis). For instance, each linear state ijj" belongs to the 
set Sn defined in Eq. ^ and ||x(s)||x = 0{\s\). It 
follows then, that the solution ip satisfies the estimate 

mmx<\s\\\r\\x + ois^), 

in the neighborhood of the bifurcation point (i?„,0). 
Therefore, since the set Sn is open, we get from the 
above estimate, that e for |s| < J. (b) Now for 
all (/i, ip) € Csr, and each n = 0, 1, . . . , if — 1, we con- 
sider the indicator function 



/(M,V') 



1, if ?A G Sn, 

0, if -0 = Oi f- = E„i, m n. 



that is, f{ii, 0) = if the branch Ce„ meets the axis 
(/i, 0) in another eigenvalue E„i ^ En. Note that / is well 
defined due to the two possibilities described by Theorem 
[H From (a) we have that if (/i, 0) is in some small neigh- 
borhood of {En, 0), then /(/i, 0) — 1. Thus, the function 
/ is constant (and equals to 1) in a small neighborhood 
of {En , 0) , and cannot change value in this small neigh- 
borhood, i.e., / is locally constant. The set Sn is open 
and the function / is locally constant on the connected 
set Ce„. Both facts clearly imply that / is continuous. 
Therefore, f{CE„) should be also connected, since the 
image of a connected set through a continuous function 
should be connected. However, / is integer valued, and 
the fact that f{CE„) is connected, implies that / should 
be constant, f — 1, for all (/i, ^) G Ce^^. Therefore, C^^ 
cannot contain a point {Em , 0) with Em ^ En and 
should be unbounded. 

Concerning the concavity of the branch, due to Propo- 
sition [21 each branch Ce„ is concave up at least for 
|s| < i5. To prove that is not possessing maximum or 
minimum points, we will apply a contradiction argument. 
Let us assume that the branch Ce„ has a local maxi- 
mum at some point. Then, due to the C^-property of 
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the branch C^;^, and since the branch is connected and 
unbounded, it follows that C^^ should possess a local 
minimum. However, as it is shown in Theorem |3] in the 
Appendix, such a minimum (here possibly attained at 
some /i), can exist in the case of a DNLS Eq. ([5]) consid- 
ered in the higher-dimensional lattice iJ^ , A/" > 1, with 
power nonlinearity, namely F{z) = \z\^° z, only in the 
case cr > Hence such a minimum in the case of a 
ID-lattice can only exist when cr > 2, which is excluded 
for the time-independent DNLS Eq. ([50]) with the cubic 
nonlinearity of a = 1 . o 

We have rigorously proved that a nonlinear state of 
Eq. (IT51) can be created by a continuation of its lin- 
ear state in [i,. Our analytical results can directly be 
compared to numerical ones: indeed, using a Newton- 
Raphson method, we can construct such nonlinear states 
starting from their linear counterpart. In Fig. [31 we plot 
the number of atoms N — I'AjP of the first three 
states, namely the ground state (solid line), first-excited 
state (dashed line) and second-excited state (dotted line) , 
as a function of the chemical potential [i. The corre- 
sponding branches begin from the linear limit (N = 0), 
where [i equals the energy of the pertinent linear state, 
and are concave up, in accordance to the analysis pre- 
sented above. The insets of Fig. |3] show the profiles of 
these nonlinear states for /i = 1.2. It is important to 
notice that, similarly to the continuous case^i^i^i^, the 
excited nonlinear states transform into a chain of discrete 
dark solitons: the first-excited state corresponds to a sin- 
gle dark soliton (one node in the wavefunction profile - 
see the middle inset of Fig. [3]), the second-excited state 
corresponds to a pair of dark solitons (two nodes in the 
wavefunction profile - see the right inset of Fig. and 
so on. 



B. Continuation from the anti-continuum limit 

Before discussing in detail the stability of nonlinear 
states in the form of discrete dark solitons, in this sub- 
section we will consider the existence and stability of non- 
linear states near the AC limit, in order to appreciate the 
emerging bifurcation structure. 

Near the anti-continuum limit, corresponding to lattice 
spacing a ^ cx), it is straightforward to find solutions of 
Eq. ([50)1 in the following form: 



V', =exp(z0,)^;.-if72(aj7, (48) 

where Qj denotes the phase. The density of the 

above solution resembles the density profile that can be 
obtained, in the Thomas- Fermi limit^, from the con- 
tinuum GPE, Eq. Thus, all the solutions for any 
number n of excited sites can be constructed following 
Eq. (|48p . One can find the analytical expression for the 
chemical potential with respect to the number of atoms 
for any such configuration of n excited sites: indeed, us- 




-2 2 -2 2 -2 2 -2 2 

i i j i 



FIG. 4. (Color online) Top left panel: The normalized num- 
ber of atoms N ja as a function of the chemical potential /x, 
for a = 10 (i.e., in the vicinity of the anti-continuum limit) 
and = 0.1. The black square indicates the region where 
this panel is magnified, as shown in top right panel. The let- 
ters A,B,...,H denote certain points in the diagram for which 
corresponding wavefunction profiles are shown in the middle 
and bottom panels. Stable (unstable) branches and respective 
states are depicted by solid (dashed or dotted) lines. 

ing Eq. (|33p and the solutions (|48l) , the normalized num- 
ber of atoms N ja reads: 

Nja^ny.-Y^Vi'iajf, (49) 

3 

where the sum runs over the excited sites. From the 
above result, it can easily be found that the slope 77 = 
d{N / a) / dji — n (for fixed n) does not depend on j - 
i.e., which particular sites are excited - but only on the 
number of excited sites. In the top panel of Fig. SI we 
show the dependence of N/a on the chemical potential 
^, for states consisting of up to three excited sites. Note 
that in this figure we have used the value a = 10, but 
we have checked that qualitatively similar results can be 
obtained for larger values of the lattice spacing. As it is 
observed in the figure, N/a depends linearly on ^ near 
the AC limit, in agreement with the analytical prediction 
of Eq. Notice that the latter is, strictly speaking, 

valid in the limit of a — !■ 00, but the linear dependence of 
N/a on /i persists for the chosen finite value of a, except 
at particular slope-changing critical points explained be- 
low. 

Let us now describe the result of Fig. |4]in more detail. 
We start with the (blue solid line) branch, correspond- 
ing to the simplest possible configuration, with only the 
center site (at j = 0) excited; an example of a state of 
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-0.15 



Parity-symmetric F 



0.57 
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0.59 



0.6 



N/a 



FIG. 5. (Color online) The pitchfork bifurcation, relevant to 
the (green) branches D, F (and its parity-symmetric one) and 
G, as viewed by the difference ^pi — of the two outer sites 
as a function of the normalized atom number N/a. 



this branch is shown in the first of middle panels of Fig. |3] 
(state A, with "A" in the top left panel marking the re- 
spective values of N/a and /i; a similar notation is used 
for the other branches below). This branch starts from 
the origin, with slope 77 = 1, but for values of chemical 
potential /i > /ii^^ = ^n^{aj)^ = 0.5 (for j ^ ±1) it 
changes slope, namely 77 = 3, as two more sites are ex- 
cited; an example of such a state belonging in this branch 
for fi > ^ci is state B (second middle panel of Fig. H]). 
Notice that further increase of /i results in a similar be- 
havior for higher values of j (not shown), i.e., this branch 
changes slope at characteristic values of the chemical po- 
tential /ii™'' = iil^(aj')^ (for j = ±to), as more sites are 
excited. This (ground state) branch of solutions is found 
to be stable throughout its continuation. 

Another branch, of slope 77 = 3, but with the three 
excited sites featuring an asymmetric configuration (see 
state C in the third middle panel of Fig. U) , also exists 
for values of N/a smaller than the ones pertaining to 
(blue) branch B. The states of this branch, which is de- 
picted by a dashed purple curve, are unstable with the 
corresponding excitation spectra being characterized by 
a pair of imaginary eigenvalues. The branch C coexists 
with another one, with two excited sites, namely branch 
E (depicted by a solid purple curve), which has a slope 
1] = 2. The states belonging to branch E (see, e.g., the ex- 
ample in the bottom left panel of Fig. 2]) are stable. Both 
branches C and E continue (as N and /i are decreased) 
up to a certain value of the chemical potential, /i ~ 0.555, 
where they collide and annihilate through a saddle-center 
bifurcation. Notice that still another branch of slope 
Tj = 2 exists, namely the (yellow) branch H, which starts 
from the linear limit (at = 0.5) and remains stable at 
least up to /i = 1.2, as shown in Fig. |4l The states of 
this branch are anti-symmetric - see the example in the 
bottom right panel of Fig. 

We now focus on another branch, of slope i] — 1 (with 
the sites j = excited), which exists for values of N/a 
greater than the ones pertaining to the (blue) branch 



A; an example of a state belonging to this (solid green) 
branch is state G - see third bottom panel of Fig. 21 As /i 
is decreased, the states belonging to this branch are sta- 
ble down to the value of chemical potential /i = 0.557: at 
this point, the slope rj changes sign, i.e., 77 < 0. The 
change of slope manifests instability according to the 
slope criterion, as suggested by the general stability cri- 
teria summarized in Section III of Refi^. It should be 
noted that instability occurs when either the slope cri- 
terion (well known also as a Vakhitov-Kolokolov (VK) 
criterion^), or the spectral condition~^° fails. It is in- 
teresting to remark that the Sturm-Liouville-type analy- 
sis discussed in Sections III and IV of the present paper, 
implies that the abstract set-up ^^'^° for the implementa- 
tion of the spectral conditions discussed therein, is valid 
for our problem. Another interesting observation is that 
the states of branch G are positive. The positivity sug- 
gests the validity of the abstract stability criteria"*^ for 
"positive solitons" which are applicable in NLS-type sys- 
tems with linear and nonlinear spatially dependent po- 
tentials, and are associated with VK-type slope criteria. 

Note that the excitation spectra of the states belong- 
ing to the continuation of branch G for < 0.557 are 
characterized by a pair of imaginary eigenfrequencies. 

Next, a decrease of fi (and increase of N/a) up to 
the point fi ~ 0.55 results in a pitchfork bifurcation, 
although this is less transparent in the variables illus- 
trated in the bifurcation diagram of the top right of 
Fig- m (see also below). The three branches resulting from 
this symmetry-breaking bifurcation are the asymmetric 
branch F (dashed green line), its parity-symmetric one 
- which has the same atom number N/a - and branch 
D (dotted green line). The symmetry-broken branch F 
inherits the stability of the original branch from which 
it stemmed (i.e., branch G), while the symmetric con- 
tinuation of branch G, namely branch D is, in fact, fur- 
ther destabilized, with the excitation spectra of the perti- 
nent states being characterized by two pairs of imaginary 
eigenfrequencies. The above mentioned pitchfork bifur- 
cation is clearly illustrated in Fig.[Sl where the difference 
t/ji — '0-1 of the two outer sites is plotted as a function 
of the normalized atom number N/a (the notation, in 
terms of the use of dashed and dotted lines, is the same 
to the one used in the top right panel of Fig. S]). 

It should be remarked that in general the modes dis- 
cussed in Fig. m cannot be expressed analytically. On the 
one hand, analytical expressions could be derived under 
the assumption that these structures are, in fact, com- 
pactly supported, i.e., that the solution is only supported 
on a few (e.g. three) sites. It is not hard to see that 
this is not true, by considering the equation of the "first 
vanishing" site. Hence, such an assumption is not self- 
consistent. Even if we bypass the above nontrivial con- 
cern and we assume three nontrivial sites and symmetry, 
we may inherit two cubic equations for the stationary so- 
lution elements which will result ultimately in a 6th order 
algebraic equation. Even if such an equation is solvable, 
the analytical expressions involved are so tortuous that 
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FIG. 6. (Color online) The linear stability analysis for the first-excited state, corresponding to a single discrete dark soliton. 
The three top panels show the real part of the lowest-order eigenvalues and the two left bottom panels show the maximum 
of the imaginary part of the eigenvalues, both as functions of the chemical potential ^. The bottom right panel shows the 
dependence of the critical value of the chemical potential, /ic, for the onset of the instability as a function of a. Branches 
shown with circles (in red) in the two top left panels denote dynamically unstable modes, which have emerged upon collision 
of modes with opposite Krein sign. The parameter values are a = 1.2 (left column), a = 0.6 (middle column), and a = 0.1 
(right column); the trap strength is in all cases 57 = 0.1. 



there is no significant intuition to be gained from this 
process. 

Concluding this section, it is important to notice that 
the study of the rich bifurcation structure presented 
above highlights the existence of nonlinear states (such 
as the ones corresponding to the branches C, F and E in 
Fig. 2]) without a linear counterpart. 

V. STABILITY OF THE NONLINEAR STATES 

A. The first-excited state 

As previously discussed, the ground-state of the sys- 
tem has been found in the discrete case to be always 
stable (i.e., for every value of a). On the other hand, as 
concerns the stability of the excited nonlinear states (per- 
taining to discrete dark multi-solitons as the nonlinearity 
increases), we note the following. Since we are interested 
in investigating the effect of discreteness on the stability 
of these states, we have performed the BdG analysis for 
three different values of the lattice spacing a (and fixed 
value of the trap strength, = 0.1). In particular, we 
have considered the following cases: a — 1.2 (correspond- 
ing to a strongly discrete case), a = 0.6 (corresponding 
to a moderate discreteness), and a = 0.1 (correspond- 
ing to a nearly- continuum setting); respective results are 
shown in Figs. IH] and H] for the first- and second-excited 
state, respectively. 

In the top left panel of Fig. [S] we show the real part 
of the lowest (four) eigenvalues, while in the bottom left 
panel we show the maximum of the imaginary part of 
the eigenvalues, for a — 1.2 (strongly discrete case). The 
branch indicated with circles (in red) has emerged from 
the collision of two modes with opposite Krein signs, 
namely the anomalous mode and the lowest positive en- 
ergy mode, and it is unstable up to the value = 0.5 
of the chemical potential. This unstable branch collides 



with the next mode producing no instability, but from 
the value — 0.55 onwards the anomalous mode starts 
colliding with higher-order modes, thus producing a new 
branch which is unstable for all values of /i > 0.55. Ac- 
cordingly, a small stability window is shown to form in 
the bottom left panel of Fig. [H] (for 0.5 < /i < 0.55). For 
a smaller value of the lattice spacing, a — 0.6 (i.e., for 
moderate discreteness) , the collision between the anoma- 
lous and the first positive energy mode occurs for a larger 
value of chemical potential, i.e., for /i = 0.86; thus, the 
configuration is initially stable, then it is characterized 
by an instability window for 0.86 < fi < 1.15, and it 
becomes again stable for 1.15 < /i < 1.35 (see bottom 
middle panel of Fig. [5]) . After a small stability window 
(for 1.55 < fi < 1.6), the anomalous mode continuously 
collides with higher-order modes resulting to instability 
for all values of /i > 1.6. Notice that the existence of such 
(in)stability windows was reported in Ref.^^, where a sim- 
ilar BdG analysis was performed (solely for the single 
discrete dark soliton configuration in a harmonic trap) . 

It is important to note that for an even smaller value of 
the lattice spacing, i.e., for a = 0.1 (close to a nearly con- 
tinuum configuration), the anomalous mode never col- 
lides with the first positive energy mode and, thus, the 
first-excited state (corresponding to a quasi-continuum 
single dark soliton) is always stable - see top right panel 
of Fig. [B]- in accordance with the findings of Refsi^i^. 
The same result can also be concluded by the bottom 
right panel of Fig. |6l where the critical value of the 
chemical potential for the onset of instability (i.e., for the 
collision between the anomalous and Kohn modes) is a 
monotonically decreasing function of the lattice spacing 
a: this indicates that (instability) stability is expected 
in the (discrete) continuous limit of the model. 

We conclude the study of the stability of the first ex- 
cited state by investigating the dynamics of unstable con- 
figurations. In particular, in Fig. [71 we show the evolution 
of an unstable discrete dark soliton, corresponding to pa- 
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FIG. 7. (Color online) The top panel shows a spatio-temporal 
contour plot of the density of the first-excited state (corre- 
sponding to a single discrete dark soliton) , for parameter val- 
ues fi — 1, a — 1.2, and SI — 0.1. The soliton stays at rest, 
up to to t ~ 1500, and then starts to perform oscillations 
of growing amplitude. The bottom panel shows the initial 
density profile. 



rameter values fi — I, a — 1.2, and Q = 0.1, as obtained 
by direct numerical integration of the DNLS Eq. ([5]). The 
initial condition, chosen in an unstable region with a rel- 
atively high instability growth rate (see bottom left panel 
of Fig.[6|), is a discrete dark soliton shown in the bottom 
panel of Fig. [T] As shown in the top panel of Fig. [71 the 
discrete dark soliton is at rest up to i « 1500; then, the 
instability sets in (due to the numerically-induced noise 
generation) and the soliton starts to perform oscillations 
of growing amplitude - a typical scenario occurring when 
a dark soliton is subject to an oscillatory instability (see, 
e.g., Ref4^). 



B. The second-excited state 

We proceed with the stability analysis of the second- 
excited state, corresponding to a dark soliton pair; our 
basic results are presented in Fig.|Hl First, we note that a 
fundamental difference of this case with the previous one 
is the existence of a second anomalous mode (recall that 
the number of anomalous modes in the excitation spec- 
trum equals the number of dark solitons22iiSii^). The 
first anomalous mode (the one with the smaller eigenfre- 
quency corresponding to the in-phase motion of the two 
dark solitons^^. - see the red branches in the top panels 
of Fig. [5]) follows a behavior similar to the one found in 
the single-dark soliton state. Thus, the discreteness in- 
duced instability presented in the previous case persists 
also in the two-soliton configuration. As concerns the 
behavior of the second anomalous mode (the one with 
the larger eigenfrequency corresponding to the out-of- 
phase motion of the two dark solitons'^^^ see the green 
branches in the top panels of Fig. [S]) we note the follow- 
ing. Starting with the left column panel (for a = 1.2), 



it is observed that the second anomalous mode initially 
resonates with the second positive energy mode, and an 
unstable quartet of eigenfrequencies (depicted in green) 
emerges. This quartet persists up to /i = 0.6, where the 
two modes split. This way, a stability window is formed 
(see the bottom left panel of Fig. [5]) which, however, is 
effectively reduced by the instability induced from the 
first anomalous mode; in fact the stability window corre- 
sponds to 0.7 < ^ < 0.74, an interval defined by the un- 
stable branch corresponding to the first anomalous mode 
(compare the red and green lines in the bottom left panel 
of Fig. [5]). Next, the second anomalous mode collides 
with a higher-order mode producing no instability but, 
eventually, further collisions with higher modes lead to 
instability. 

For a smaller value of the lattice spacing (a = 0.6), 
and contrary to the previously examined - highly dis- 
crete - case of a = 1.2, the first anomalous mode is 
initially stable, but becomes unstable for ^ = 0.85. On 
the other hand, the quartet that has emerged from the 
second anomalous mode and the second positive energy 
mode (which is initially unstable as before) splits at 
fi = 0.8; this way, a small stability window is created for 
0.8 < /i < 0.85 (see bottom middle panel of Fig. [5]), while 
for ^ > 0.85 the configuration is dynamically unstable. 
Finally, for a = 0.1 (corresponding to a quasi-continuum 
configuration), the right column panels of Fig. [5] sug- 
gest that an instability induced by the second anomalous 
mode occurs for fi < 0.7, but then, for fi > 0.7 the con- 
figuration remains stable (although for sufficiently large 
fi it will become unstable again). This result is in accor- 
dance with the findings of Rcfs.^^'^^, which suggest that 
in the continuum limit the multi-soliton solution is un- 
stable near the linear limit (due to the second anomalous 
mode) and may only be unstable thereafter in parametric 
windows due to collisions of the second anomalous mode 
with higher positive energy ones. 

In the left set of panels of Fig. ^ we show the dynam- 
ics of an unstable two-dark soliton configuration, corre- 
sponding to parameter values ^ = 0.5, a = 0.6, and 
il = 0.1, as obtained by direct numerical integration of 
the DNLS of Eq. The initial condition is a discrete 
two-dark soliton state shown in the bottom left panel of 
Fig. O As shown in the top left panel of Fig. [31 the dis- 
crete dark soliton pair is stationary up to i « 500; then, 
the instability manifests itself: the out-of-phase motion 
of the two-soliton state excites the second positive energy 
mode (often referred to as quadrupole mode in the con- 
tinuum case-) of the system. This excitation results in 
a breathing behavior of the configuration. This type of 
instability was also found in the continuum counterpart 
of the system (see Fig. 5(c) of Refi^). Note that this 
particular simulation corresponds to the first instability 
band shown in the bottom middle panel of Fig. [51 The 
results of simulations performed with values of fi corre- 
sponding to the second and third instability bands, i.e., 
for /I = 1 and fi = 1.5, are respectively shown in the mid- 
dle and right panels of Fig. [51 In both cases, the initially 
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FIG. 8. (Color online) The BdG analysis for the second-excited state, corresponding to a discrete dark soliton pair. The 
top (bottom) panels show the real (imaginary) part of the lowest-order eigenvalues as functions of the chemical potential ^l. 
Branches shown with circles (in red or green) in the top panels denote dynamically unstable modes, which have emerged upon 
collision of modes with opposite Krein sign. The parameter values are a = 1.2 (left column), a = 0.6 (middle column), and 
a = 0.1 (right column); the trap strength is in all cases SI = 0.1. 




FIG. 9. (Color online) Same as Fig. [T] but for the second-excited state, corresponding to a discrete dark soliton pair, for 
parameter values fj, — 0.5 (corresponding to the first instability band - see bottom middle panel of Fig. [8|, a = 0.6, and 

n = 0.1. 



quiescent two-dark-soliton configuration becomes unsta- 
ble and is set into motion, with the solitons performing 
an in-phase motion. It is clearly observed that the insta- 
bility manifests itself at different times in the two cases, 
namely at < « 800 and t « 180 for the middle and right 
panels, respectively; furthermore, the amplitude of os- 
cillation of the dark soliton configuration in the former 
case is much smaller than the one shown in the latter. 
Thus, although both cases correspond to linearly unsta- 
ble two-dark-soliton configurations, the one shown in the 
middle panel appears to be more robust than the one in 
the right. This may be partially connected to the insta- 
bility growth rates, but perhaps, more importantly, to 
the different modes of the background with which the 
internal in-phase soliton mode resonates in the different 
cases. 



VI. CONCLUSIONS 

In this work, we presented a systematic study of the 
existence, stability and bifurcations of nonlinear states 



of a self-interacting quantum harmonic oscillator (QHO) 
on a lattice. The considered model, namely a discrete 
NLS equation incorporating a (discrete) harmonic trap, 
may be used to describe the dynamics of an array of 
BEC droplets in a deep optical lattice, but also discrete 
nonlinear guided- wave optical systems. 

Our considerations started with the analysis of the per- 
tinent linear problem. We presented the energy spectrum 
and the eigenstates of the linear problem as functions of 
the lattice spacing a. This way, we spanned all pos- 
sible cases, starting from the continuum limit (i.e., the 
well-known QHO for a 0) to the anti-continuum one 
(a — oo), where the ground state energy asymptotes to 
zero, while the excited states exhibit a parabolic energy 
spectrum. Next, using global bifurcation theory (and em- 
ploying, in particular, a functional-analytic theorem from 
the work of Rabinowitz), we rigorously proved that - in 
the discrete regime - all eigenstates of the linear problem 
can be continued to nonlinear ones, so that each linear 
state possesses a nonlinear counterpart. Using this re- 
sult, we were able to construct numerically the nonlinear 
states emerging from their linear siblings; this way, we 
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found the ground state of the system (which acquires the 
Thomas-Fermi profile in the anti-continuum hmit), while 
the excited states take the form of a chain of station- 
ary discrete dark solitons. The ant i- continuum limit was 
studied as well; it was found that the solutions present 
a complex bifurcation structure, which was elucidated 
along with the stability of the corresponding branches. 
The pertinent bifurcation diagram also revealed the ex- 
istence of nonlinear states with no linear counterpart. 

We also performed a detailed linear stability analysis of 
the different nonlinear solutions ensuing for different val- 
ues of the lattice parameter a, solving the BdG equations 
eigenvalue problem, for the ground state, as well as for 
the first and second excited states (the latter, correspond 
to a single dark soliton and a pair of dark solitons, re- 
spectively, in the strongly nonlinear regime). While the 
ground state was found to be completely stable for all 
values of the lattice spacing, this was not the case for the 
discrete dark soliton states, which revealed a quite rich 
stability spectrum. In the strongly discrete regime, the 
single dark soliton was found to be potentially unstable, 
due to collisions of the first anomalous mode eigenvalue 
with the rest of the normal modes of the system, for in- 
creasing chemical potential ^. As the system becomes 
more continuous, i.e., for decreasing lattice spacing a, 
stable windows appear and gradually expand; eventually, 
in the quasi-continuum regime, the anomalous mode re- 
mains below the positive energy mode for all values of /x, 
and the soliton becomes stable. A similar behavior (from 
the strongly discrete to the quasi-continuum regime) but 
with additional sources of potential instabilities (from 
the additional anomalous mode) was identified for the 
two dark soliton configuration. The anomalous mode re- 
sponsible for out-of-phase motion between the solitary 
waves is initially in resonance with the second positive 
energy mode, thus creating an instability, but eventually 
they split to create small windows of stability (in the 
discrete regime). In the quasi-continuum limit, the first 
anomalous modes yields no instabilities while the second 
one leads to windows of instability (which are more pro- 
nounced near the linear limit). 

In both cases, our analysis revealed that the effect of 
discreteness is to chiefiy offer additional sources of insta- 
bility of the dark soliton states for atom number param- 
eter ranges for which they would be in the continuum 
counterpart of the model. This is due to the pronounced 
dependence of soliton anomalous modes on the chemical 
potential, as well as due to the discreteness eliminating 
some of the symmetries (such as the dipolar symmetry of 
the first positive energy mode) present in the continuum 
limit. 

We would like to conclude by mentioning some main 
differences between the results concerning ([T^ and ([50]). 
and the DNLS equation without potential^. One of 
these differences concerns the infinite lattice limit, es- 
pecially in terms of the bifurcation analysis carried out 
herein. The case of the parabolic potential retains its 
point-spectrum nature even in that limit and the spac- 



ing of the energy levels is chiefly controlled by the trap 
frequency Jl. On the contrary, in the absence of 
as the lattice becomes infinite in the realm of the above 
paper, the point spectrum due to the finiteness of the 
domain converts itself into a continuous spectral band 
and hence its properties (and bifurcations) are substan- 
tially different. The bifurcation mechanism analyzed 
herein (bifurcation from simple eigenvalues) is one of the 
main types for generation of nonlinear states, as it has 
been highlighted in Ref.5? (see Section 11^, pg. 046602- 
2). Another relevant difference concerns the lengthscales 
of these states. In the absence of fl, the point spec- 
trum eigenfunctions are spatially " extended" (within the 
length-scale of the lattice). On the other hand, the spa- 
tial eigenfunctions in the case of the parabolic trap prob- 
lem are localized within a lengthscale controlled by the 
trap frequency fl. Hence, the presence of a parabolic trap 
yields an additional lengthscale which can be used to in- 
duce interesting phenomena, such as for example the ones 
that emerge from the competition of the trap lengthscale 
with the intrinsic lengthscale of the lattice. This is e.g. 
what produces the complex bifurcation diagrams such as 
the one of Fig. SI while such a phenomenology is likely 
more limited in the context of a 2- or a 3-site lattice 
(only). 

It would be interesting to extend our considerations 
in other settings, such as ones involving different types 
of trapping potentials, multi-dimensional one-component 
systems (e.g., in the case of both dark soliton and vortex 
type entities in two-dimensional settings), as well as in 
multi-component systems. Work is in progress in these 
directions and relevant results will be presented in future 
publications. 
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APPENDIX 

In this section we give a proof on the existence of an 
excitation threshold in the sense of Refi^, for the DNLS 
equation ^ considered in the lattice Z^, Af > I. For 
technical purposes it is more convenient to work with the 
focusing version of DNLS of Eq. ([5]) , having the opposite 
sign on the nonlinearity. We shall reduce the DNLS of 
([S]) to the one with an effectively opposite coefficient of 
the nonlinearity under the, so-called, staggering trans- 
formation. This transformation is defined as (see, e.g.. 
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Ref, 



50^ 



(50) 



for j :— {ji, j2, ■ ■ ■ , jj^) G ^■'^ (a trivial multiplication 
by a suitable phase factor is also needed to form the 
corresponding local term within the discrete Laplacian). 
Thus, under Eq. (j50p . the A/'-dimensional, focusing ver- 
sion of the system ([5]) with a general power nonlinearity, 
can be actually written as (taking advantage of 
the time-reversal symmetry) 

i^, -f ^A2^, + ia^rf^ljlV, -f |V,f ^V^, = 0, (51) 

considered in the A/'-dimensional cube of l/^ with edges 
of length 2L, 

2 = {{xn^-'-^^iM) ■ < Ji,---,jV <K + l}, 
2L 

Xj,^-L+j,a, a= 1 = 1,..., A/. 

K + L 

The discrete eigenfunctions on Q are denoted by 
The interior of the cube Q is given by 



) : l<Ji,---,M<K}, 



and (|51|) is supplemented with Dirichlet boundary con- 
ditions 



ipj=0, ondQ := Q\ Q. 



(52) 



Solutions ipj — )■ ipj exp{—ifit), of (|5ip -([52 |) . are equiva- 
lently, solutions of the constrained minimization problem 



In - {H[ij] : N[ij] = 7^} 



(53) 



where the chemical potential /i appears as a Lagrange 
multiplier associated to the minimizer of (|53p . In ()53p . 
denotes the Hamiltonian 

urn = ^(-^2f,i')2 + lc'a'Y.ijfi*ii' 



2a+2 
3 1 ' 



while 



Q 

the norm of the Hilbert space P of square summable 
sequences, represents the atom number or optical power 
(see ([55)) ). We have the following 



Proposition 3 I-ji > i/ and only if TZ satisfies the 
inequality 



, (54) 



for all tj} G 



K+2) 



Proof: By the definition of I-ji in ([55]) it follows that 
Ifl > if and only if 

s 

+ ia^l7^X]|jn^,P, (55) 

e 

for aU V e M-^(^+2). Let now V e M-^(^+2), ^ ^ ar- 
bitrary, and consider the element z — ^/iZW^W^'^'tp. Ob- 
serving that N[z] = ||z||| = R, by substitution of z in 
(155]) we derive dMl). o 

From the Proposition 4.2, p. 680 of Ref.^^, it clearly 
follows that if (T > , there exist a constant C > such 
that for any e > 0, the inequality 



J62 



x^(-A2V,V')2, (56) 

holds for aU ipee'^. Thus, it is an immediate consequence 
that there exist C > 0, such that 



E m'^^'<c\ E i^.i 



(57) 



for all ij: € P and cr > We define for brevity the 
functional 

^= i(-A2^, V)2 + ia'f^' E l^fl^^-l'- 



In analogy with Eq. (4.2), p. 680 of Ref.^, if C* is the 
infimum over all the constants for which ([57]) holds, then 
C* it can be characterized as 



inf 



(58) 
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Therefore, the excitation threshold 7?.throsh for the DNLS 
equation i51\) wih be defined by a comparison of (j54p and 
(1571) as 



(59) 



(cr + l)(72.thrcsh) = C*. 

We conclude with tlic following theorem. 



Theorem 3 Let ^ > jj 

11 = 



A. Assume that 



TZ. Then 



1 - 



n 



thresh 



(60) 



B. If TZ < 7?.thrcsh then T-jz = and there is no ground 
state minimizer of 153)) . 

C . If TZ > 7?.thresh then T-ji < and there exists a mini- 
mizer of the variational problem ^5S^] . 



Proof: A. Let us note first that ((56)) holds for any 
-0 e with the same optimal constant C*, since 

^N'{K+2) jg ^ finite dimensional subspace of P. Then, 
using (j56p with its best constant C» we derive that 



1 ^ 

Q 



> E[,p] - (7^thresh)""7^-^^[V'], 

thus ([eo]) . 

B. Assuming that TZ < 7?.thrcsh, it follows from (|60l) that 
I^TZ ^ 0- On the other hand, we may consider some ip G 



where A > arbitrary. 



7^ 



A 



Considering the element zx — V^||V'||2 ^"0 we observe 
that 



and 



\zx\\l=TZ 



H[zx]^X'E[^]~^Y.\'f'^ 

Q 



2a+2 



For A sufficiently large, we get that II[z\] < 0. Therefore 
if 7^ < T^thresh we should have 2-jz = 0. Assuming that 
this infimum is attained at a state 0, then X-ji = implies 
that 



E[4>]^ 



1 



\4>n 



\2a+2 



Q 

N[$] = Y,\^,\'=TZ. 
Q 



(61) 



Then, inequality (j56l) with its best constant C,, if in- 
serted into (|6ip . is giving the contradiction 



El 



< 



1 



I 2 |2(T+2 - 1 



1 - 2a2 



7^ 



thresh 



<i?[<^]. 



C. By the definitions dSS]), (|59l) of C* and 7^throsh re- 
spectively, it follows that if 7?, > 7?.thicsh then a 0* € 
should exist which does not satisfy inequality ([54]) . hence 
I-jz < 0. Indeed, such a minimizer exists since Hlip] is 
bounded from below and in the finite dimensional space 
^Ar(K+2) ^Yie infimum I-r < is attained, o 

Let us note that the complementary results presented 
in this section are of independent interest since they prove 
existence of nonlinear states for the DNLS Eq. ([5|) di- 
rectly, together with the existence of a threshold for their 
activation energy. It is important to note that in the con- 
tinuous limit a — ^ oo, an excitation threshold exists only 
in the critical case cr = jj (see Sections 3 and 4 of Refi^). 
We also remark that the results can be extended in the 
case of the infinite lattice l/^ by implementing the con- 
centration compactness arguments^^^. 
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